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ABSTRACT 

We perform a large set of cosmological simulations of early structure formation and 
follow the formation and evolution of 1540 star-forming gas clouds to derive the mass 
distribution of primordial stars. The star formation in our cosmological simulations 
is characterized by two distinct populations, the so-called Population III. 1 stars and 
primordial stars formed under the influence of far-ultraviolet (FUV) radiation (Popu¬ 
lation III.2 d stars). In this work, we determine the stellar masses by using the depen¬ 
dences on the physical properties of star-forming cloud and/or the external photodis- 
sociating intensity from nearby primordial stars, which are derived from the results 
of 2D radiation hydrodynamic simulations of protostellar feedback. The characteris¬ 
tic mass of the Pop III stars is found to be a few hundred solar masses at z ~ 25, 
and it gradually shifts to lower masses with decreasing redshift. At high redshifts 
z > 20, about half of the star-forming gas clouds are exposed to intense FUV radia¬ 
tion and thus give birth to massive Pop III.2 d stars. However, the local FUV radiation 
by nearby Pop III stars becomes weaker at lower redshifts, when typical Pop III stars 
have smaller masses and the mean physical separation between the stars becomes large 
owing to cosmic expansion. Therefore, at z < 20, a large fraction of the primordial gas 
clouds host Pop III. 1 stars. At z < 15, the Pop III. 1 stars are formed in relatively cool 
gas clouds due to efficient radiative cooling by H 2 and HD molecules; such stars have 
masses of a few xlO M 0 . Since the stellar evolution and the final fate are determined 
by the stellar mass, Pop III stars formed at different epochs play different roles in the 
early Universe. 

Key words: methods: numerical - stars: formation - stars: luminosity function, mass 
function - stars: Population III - dark ages, reionization, first stars. 


1 INTRODUCTION 

Primordial stars play vital roles in the early evolu¬ 
tion of the cosmos by initiating cosmic reionization and 
chemical enrichment of the inter-galactic medium (e.g., 
iBromm fe Yoshidali201ll l. Stellar evolution and death, which 
regulate the dynamical, radiative, and chemical feedback to 
the surrounding medium, are largely determined by the stel¬ 
lar mass. Knowledge of the characteristic mass of primordial 
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stars, or their initial mass function (IMF), is thus essen¬ 
tial for understanding the initial stages of the first galax¬ 
ies’ formation. Whereas direct observations can be utilized 
to derive the stellar IMF in the present-day Universe, only 
limited and indirect observational constraints are available 
regarding the characteristic mass of the first generation of 
stars. For instance, the elemental abundance patterns of 
Galactic extremely metal-poor stars provide signatures of 
nucleosynthesis in primordial stars, from which the mass of 
the progenitor stars can be inferred (e.g., Cafiau_etaL 201fJ; 
lAoki et al.ll20I4l ; iKeller et al.ll20I4l ; iTominaga et al. 2014h . 

Theoretical studies offer a possibility to reveal the na- 
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ture of primordial sta rs. The well-established standard cos- 
mological model (e.g.. IPlanck Collaboration et al.ll2014l ) al¬ 
lows us to take an ab initio appr o ach to study primordial 
star formation (see B rom ml 120131: [g] ove 3 20131. for recent 
reviews). According to these recent studies, primordial star¬ 
forming gas clouds are formed in the so-called mini-haloes 
with masses of 10 5 — 10 6 M@ at very early epochs. At the 
centre of such a cosmological mini-halo, the gas cloud gravi¬ 
tationally collapses via H 2 and/or HD molecular cooling, so 
that a tiny (~ 0.01 Mq) embryo - protostellar core - forms 
(e.g.. lYoshida et al.ll2008h . 

The newly-born protostellar core grows in mass via 
rapid gas accretion from the surrounding envelope. The typ¬ 
ical mass of the accretion envelope is very close to that of the 
original gravitationally-unstable gas cloud (so-called Jeans 
mass), 

“ 1000 (i55Kr (li ’ (1) 

where T is the gas temperature and tin is the hydrogen 
number density (lAbel et al.l 12002 '). The final stellar mass 
critically depends on how long the gas accretion continues. 
A key process in this later accretion stage is the protostel¬ 
lar radiative feedback, which se ts the final ste l lar m ass by 
terminating the mass accretion. iMc.Kee fe Tanl d200Sl J eval¬ 
uate the potential impacts of the ultraviolet (UV) radiative 
feedback using a semi-analytical model. They consider that, 
when a protostar becomes massive enough to emit a copious 
amount of ionizing photons, an H n region grows in polar di¬ 
rections of a circumstellar disc. The stellar ionizing photons 
heat up and photoevaporate the gas on the disc surface. 
More recently, the feedback process has been studied with 
multidimensi onal radiation hydrodynamic (RHP) numerical 
simu lations dHosokawa et al.l 1201 ll : IStacv et al.l 120121 : ISusal 
l2013lj . demonstrating that UV radiative feedback does in¬ 
deed shut off the mass accretion and the resulting stellar 
masses can be a few tens of solar masses. 

These previous studies calculate only several individ¬ 
ual cases and thus it remains unclear if they represent fidu¬ 
cial cases. Obviously, statistical studies are needed to study 
the overall mass distribution o f primordial stars. Our pre¬ 
vious work (jHirano et alJlioiil . hereafter Paper I) achieves 
this goal with more than one hundred cosmological sam¬ 
ples of primordial star-forming clouds. We first follow the 
early evolution until the formation of the protostellar core 
with three-dimensional (3D) cosmological simulations. We 
then study the subsequent evolution in the mass accretion 
stage under the influence of protostel lar UV feedbac k by per¬ 
forming local 2D RHD simulations dHosokawa et al.ll201lf) . 
The protostellar evolution is calculated self-consistently by 
num erically solving the int erior structure of the protostar 
(e.g.. lOmukai fe Pallall2003l ). A wide range of stellar masses 
is obtained, extending from 10 to 1600 Mq , which suggests a 
great diversity of primordial stars. Despite several physical 
effects that are not directly realized in 2D simulations, e.g., 
disc fragmentation and stellar multiplicity, such diversity is 
consistently reported bv othe r studies including recent 3D 
simulations (ISusa et al.ll20lj '). Q Interestingly, we find cor- 


1 We adopt the so-called alpha-viscosity whose value is cali¬ 
brated with respect to the results of recent 3D simulations (see 


relations between the final stellar masses and the physical 
properties of star-forming clouds, 
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where Mj ea ns is the gas infall rate at the Jeans scale, as 
well as between the final stellar masses and the physical 
properties of the mini-haloes, 
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where 2 is the forming redshift and M v n is the virial mass 
(these equations are given in eqs. 13 and 19 in Paper I). Such 
correlations are useful because one can estimate the final 
stellar mass using the early physical state of the clouds or 
of the haloes, without following the detailed evolution from 
the cloud collapse up to the termination of mass accretion 
on to the protostellar core. 

The ultimate goal of our statistical study is to directly 
calculate the primordial IMF theoretically. Clearly, a num¬ 
ber of improvements over previous studies are needed to 
achieve this. First, in Paper I, we select star-forming gas 
clouds from a few different simulations (realizations) with 
small cosmological volumes. It would be desirable to use a 
single large volume to make a cosmologically representative 
sample. Secondly, it is often assumed that the primordial 
stars are unaffected by any external feedback from other 
nearby stars (Population III. 1 stars @). However, primordial 
stars forming under the influence of a UV radiation back¬ 
ground, the so-called Population III.2 stars, are thought to 
have different characteristics from Pop III. 1 stars, because 
the abundance of H 2 and HD molecules, which determine 
the thermal evolution during the cloud collapse, are vulner¬ 
able to the external radiation. The relative occurrence of the 
two populations is not known. 

In this paper, we consider a comoving volume of (3 
/i -1 Mpc) 3 for our cosmological simulations, rather than the 
simulation volume of Paper I. Furthermore, we classify the 
Pop III.2 stars into two subclasses depending on the hard¬ 
ness of the radiation: (1) the photodissociation-dominated 
case (Pop III. 2d) and (2) the photoionization-dominated 
case (Pop III.2i), where the subscripts D and I stand for 
dissociation and ionization, respectively. In the former case, 
because H 2 and HD molecules are destroyed by photodisso- 
ciating photons, the temperature in a coll apsing cloud is 
higher than for the Pop III. 1 case (e.g., Omukai fe Pallal 
l200ll ; lOmukai fe Yoshiill2003l : lO’Shea fe Normanll2008ll . The 
mass accretion rate on to the protostar increases follow¬ 
ing the well-known scaling relation M oc T 3 ^ 2 . At the 
higher accretion rates, the protostar has a larger radius and 
lower effective temperature for a given stellar mass (the 
“bloating” phenomen on, see e.g. IZinnecker fe Yorka 120071 : 
iHosokawa et alj l201ll l. This weakens the stellar UV feed¬ 
back, and the resulting final stellar mass becomes higher 
than in the Pop III. 1 case. 


appendix A in Paper I). Thus our 2D calculations include effi¬ 
cient angular momentum transfer via the non-axisymmetric disc 
structure. 

2 In Paper I, we have classified Pop III. 1 stars into two sub-classes 
depending on whether HD molecular cooling affects the thermal 
evolution of the cloud (Pop III.Ihd) or not (Pop III.1 h 2 )- 
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In the latter case (Pop III.2i), radiation fields affect the 
thermal evolution in the opposite direction; photoionization 
first enhances the abundance of free electrons, which are the 
catalyst for generating H 2 molecules, and then a collaps¬ 
ing cloud evolves at low temperature due to the enhanced 
H 2 and HD formation and cooling ('e.g.. lYoshida et al.ll2007l ; 
iHosokawa et al.l2012l 'l . Overall, these theoretical studies sug¬ 
gest that the mass of a primordial star depends critically on 
the local strength of the radiation emitted by other stars. 

Clearly, we need to consider all sub-populations in or¬ 
der to calculate the primordial IMF. In the present study, 
we consider both Pop III. 1 and III. 2d stars but ignore the 
contribution of Pop III.2i stars (see Sec. [7] for a discussion 
of this topic). In order to generate non-biased cosmological 
samples of star-forming clouds, we first choose 1540 haloes in 
our large cosmological simulation. We then derive the corre¬ 
lations between the stellar mass and the physical properties 
of the parent gas cloud to estimate the mass of a star form¬ 
ing in each cloud. To this end, we investigate how the final 
stellar mass depends on the strength of the photodissociat- 
ing radiation field by performing computationally intensive 
RHD simulations for 45 representative cases. We then ap¬ 
ply the derived formulae for the final stellar mass to our 
cosmological samples. For each cloud, the local intensity of 
the photodissociating radiation is calculated by summing 
up contributions from nearby primordial stars distributed 
in the cosmological volume. This approach provides us with 
mass distributions of Pop III. 1 and III.2 d stars at different 
epochs. 

The remainder of the paper is organized as follows. In 
Section 0 we describe the numerical methods employed in 
our set of simulations: cosmological simulations, RHD sim¬ 
ulations, and the post-processing calculations of the local 
far-ultraviolet (FUV) intensity. In Section [3] we first sum¬ 
marize our main results. Section [4] shows the statistics of 
our cosmological samples of the primordial star-forming re¬ 
gions. Section [o] presents the results of the local RHD sim¬ 
ulations, with which we also explain how the stellar mass 
is assigned to each primordial cloud. Section [6] shows the 
resulting mass distributions of primordial stars, classifying 
our statistical samples into Pop III. 1 and III. 2d cases. We 
finally give concluding remarks in Section [7] 

Throughout this paper, we adopt the standard A-cold 
dark matter (A-CDM) cosmology with the total matter den¬ 
sity ff m = 0.3086, the baryonic density fib = 0.04825, 
the dark energy density IIa = 0.6914 in units of the crit¬ 
ical density, the Hubble const ant h = 0. 6777, and the 
prim ordial index n s = 0.9611 dPlanck Collaboration et al.l 
2014|). The power spec tra are defined by the equation in 
Eisenstein fc Hul (119990 and normalized to <T8 = 0.8288. 


2 METHODS 

2.1 Sampling Star-Forming Regions in 
Cosmological Simulations 

First, we perform a large parent cosmological simulation 
with relatively low spatial resolution and identify dark mat¬ 
ter (DM) haloes. We then recalculate the evolution using 
the so-called hierarchical zoom-in technique for 55 cradle 
regions which contain a lot of haloes. From these zoom-in 


simulations, we pick out the star-forming clouds and calcu¬ 
late the subsequent collapse until the central density reaches 
10 7 cm -3 . The cloud state at this moment is used to infer the 
final stellar mass (see Sec. 14.411 . Here, we describe each step 
of our procedure, emphasizing newly developed techniques 
and the improvements over previous studies. 

Cosmological simulations are performed using the par¬ 
allel Y-body / s moothed particl e hydrodynamics (SPH) 
solver gadget-2 (ISpringell I2005T ) in its version suitably 
adopted for primordial star formation. We solve chemical 
rate equations for 14 primordial species (e - , H, H + , H - , 
He, He+, He++, H 2 , H+, D, D+, HD, HD+, HD - ) as in 
lYoshida et alJ (120061 . l2007f). We have updat ed the cooling 
rates for H 2 and HD fsee lcall i fc Pal lal 20131') as well as the 
three-body H 2 formation rates (lForrevll2013al lbh . 

We generate the cosmological i nitial condit ions using 
the modified version of N-GENIC (ISpringell 120051 1. In order 
to achieve sufficient spatial resolution to resolve the primor¬ 
dial gas clouds with masses of ~ 1000 Mg, we use a hierar¬ 
chical zoom-in technique. First, we run the parent Y-body 
(DM-only) simulation with a large cosmological volume of 
Tbox = 3 /i -1 Mpc (comoving) on a side employing 768 3 


DM particles. The particle mass of the dark matter compo¬ 
nent is 7400 Mg and the gravitational softening length is 
set to be 115 pc. We follow the evolution from a = 99 to 
9. The background panel in Fig. [Q shows the final output 
of this simulation. Primordial clouds form within the dark 
matter haloes located at the densest parts of the cosmic 
filaments. We then pick out all candidate sites for primor¬ 
dial star formation in the simulated volume. We identify DM 
haloes by running the Friends-Of-Friends (FOF) halo finder. 
Fig. [2] shows the obtained mass function of DM haloes. Our 
simulation result agrees well with the analytical prediction 
(Press-Schechter mass function), confirming that our parent 
simulation represents a typical cosmological volume. Note 
that our halo selection should be non-biased if all of the 
star-forming sites are selected from this simulation volume. 

Fig.HJshows that the dark haloes are strongly clustered. 
We locate 55 zoom-in spherical regions with the radius of 390 
h -1 kpc to cover such cradle regions (one of them is marked 
with the red circle in Fig. [T]). We only use the inner 3/4 vol¬ 
ume of each zoom-in region to look for gas clouds in order to 
exclude non-physical haloes which appear in the outer part 
of the zoom-in regions. The net volume considered as pri¬ 
mordial star-forming regions is about 20 per cent of the total 
simulation volume. We generate 55 zoomed initial conditions 
with baryonic components (SPH particles) with increased 
mass resolution by a factor of 4 3 . The resulting effective 
number of particles in a zoom-in region is 3072 3 . The corre¬ 
sponding small-scale density fluctuations are added suitably 
as given by our adopted A-CDM cosmology. The masses of 
DM particles and of the gas particles within the zoom-in 
region are 116 and 19 Mg, respectively. The gravitational 
softening length is set to be 29 pc. We continue the zoom-in 
simulations down to redshifts 2 ~ 30 — 25, when the first 
gravitationally collapsing cloud is formed. 

When non-linear objects are formed, we further clip the 
densest part out of the zoom-in regions with the following 
cut-out criteria: 


3 Other length scales described below are all in comoving scale. 
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Figure 1. Schematic view of the cosmological simulation. From a parent simulation (the background panel), we generate 55 zoomed 
initial conditions which cover each cradle region (the red circle shows one of them). We finally obtain 1540 star-forming clouds from 
zoom-in simulations (white circles in the foreground panel) and follow their formation and evolution. The colour contour reveals the 


and dynamical feedbacks. Primordial clouds forming at the 
centres of such haloes will host Pop III. 1 or III. 2d stars. 

Finally, we restart the simulations for each selected gas 
cloud. When a cloud collapses, we continually increase the 
spatia l resolution b y adopting the particle-splitting tech¬ 
nique (IKitsionas fc Whitworthll2002l ) to insure that the local 
Jeans length is always resolved by 15 times the local smooth¬ 
ing length of the SPH particles. We follow the collapse of 
each cloud until its central density reaches 10 7 cm -3 , which 
ultimately results in 1540 primordial star-forming regions 
within one parent cosmological volume. 

In the present study, we do not consider the gas metal 
enrichment and hence the formation of Pop II stars. Previous 


projected density level of DM which increase from blue to yellow. 


1. Because we stop each zoom-in simulation when a cloud 
begins to collapse, the subsequent evolution of the other 
clouds is not followed in the same single simulation. We tag 
potential sites where other clouds are later formed, by iden¬ 
tifying the regions with densities a few times greater than 
the cosmic mean value. These zoomed-in regions normally 
have ~ 10 — 100 candidates. 

2. We select haloes which lie in the densest region within 
1 physical kpc radius around the density peak, implicitly 
assuming that the halo forming there will not be affected 
by external feedback from nearby sources except the pho- 
todissociating feedback, including photoionizing, chemical, 
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Figure 2. Mass function of DM haloes. The crosses show the re¬ 
sult of the parent simulation at z = 9. We run the FOF halo finder 
with linking length b = 0.2. We retain haloes containing more 
than 50 IV-body particles. The solid lines show the analytically 
calculated Press-Schec hter mass function at z = 9 (thick), 20, 
and 30 (thin) using the lReed et al.1 lldOOil ) code. The dashed hor¬ 
izontal line represents the number density of 1 [(3 h~ 1 Mpc) -3 ], 
with which one halo should be found in the cosmological volume 
of the parent simulation. 


studies show that the FUV radiation from nearby Pop II 
stars dominates over the Pop III component at z <10 — 15 
(e.g., lAgarwal et al.l [20121 : I Johnson et al.ll2013l ). Thus, we 
stop our simulations at 3 = 10 (see also discussion in Sec. 

m - 

2.2 Deriving the Correlations between the FUV 
Intensity and Stellar Masses: 2D RHD 
Simulations 

The mass of a Pop III. 1 star is largely determined by the 
physical properties of the natal gas clouds, e.g., the cloud 
mass and the rotation degree as shown in Paper I. For 
Pop III. 2d stars, however, there are additional parame¬ 
ters that could affect the star formation process. One is 
the intensity of photodissociating FUV radiation in Lyman- 
Werner (LW) bands, which is often normalized as J 21 = 
Jlw/(10 21 erg s _1 Hz -1 sr -1 ). This radiation destroys H2 
and HD molecules, and prevents the cooling and collapse of 
the cloud. In this study, we investigate the dependence of 
Pop III. 2d star-formation on the parameter J 21 . 

Another parameter is the central density when the pho¬ 
todissociating photons reach the cloud. We fix this param¬ 
eter at 7iH,cen,rad = 10 cm” 3 for the following reasons. In 
reality, the FUV photons might affect the earlier phase 
of the collapse (nH,cen,rad < 10 cm -3 ). Pop III. 2d stars 
could then form with less efficient self-shielding. However, 
in our calculations it takes ~100 Myr for such a low-density 
cloud to begin to collapse with the reduced amount of H 2 
molecules. It is unlikely that the intensity of FUV radiation 
remains strong enough for the Pop III. 2d star formation 
mode for such long time; typical stellar lifetimes of the FUV 
sources are only a few million years. For the opposite case 


(««,eeryrad > 10 cm -3 ), self-shielding prevents photodisso¬ 
ciation and the cloud forms a Pop III. 1 star; we find that 
self-shielding becomes sufficiently efficient in the cloud at 
densities ~ 10 2 — 10 3 cm -3 . We thus adopt the critical value 

^H,ceil, rad — 10 Cm 

We calculate the evolution of nine different clouds irra¬ 
diated by photodissociating photons. We study the effects of 
external irradiation, varying the intensity as a free param¬ 
eter with J 21 = 0, 0.1, 0.316, 1, and 10. We have explicitly 
checked that for J 21 = 0.01, the thermal evolution of the 
collapsing cloud is nearly identical to the Pop III. 1 case (J 21 
= 0). On the other hand, the collapse is almost completely 
prohibited for intensities exceeding J 21 = 10. We thus ex¬ 
amine the parameter range of J 21 = 0.1 — 10. 

We study the formation of Pop III. 2d stars in a two- 
step manner as follows. First, we calculate the evolution 
of collapsing clouds under the FUV fields by gadget-2. 
A uniform FUV radiation field is assumed here. We em¬ 
ploy the self-shielding functions (I Wolcott-Green fc Haimanl 
l201ll : [Wolcott- Green et al.ll201ll ) for H 2 and HD molecules. 
For each SPH particle, we calculate the column densities 
along six directions (±X, ±Y, ±Z) to account for the di¬ 
rectional dependence of the self-shielding effect. We stop 
the collapse simulations when the central density reaches 
n H.cen = 10 13 cm -3 and assume that an optically thick hy¬ 
drostatic object (protostar) has formed at this point. The 
final configuration is used to generate the initial conditions 
for our 2D RHD grid-based simulations of the subsequent 
protostellar accretion stage by azimuthal averaging around 
the rotational axis ((^-direction). 

We follow the subsequent evolution after the birth of 
the protostar by conducting 2D axisymmetric RHD simula¬ 
tions of the collapsing cloud coupled with the simultaneous 
stellar evolution of t he central accreting protost ar located 
in a central sink cell dHosokawa et al.| [2011. 2012, Paper I). 
We follow the evolution until the mass accretion rate falls 
below 10 -4 Mq yr _1 in each case. We regard the mass at 
this moment as the final stellar mass. Whereas we normally 
follow the accreting protostar’s evolution by solving for the 
stellar interior structure, it occasionally becomes difficult 
to construct a stellar model with highly variable mass ac¬ 
cretion histories. This often happens in Pop III. 2d cases, 
when the mass accretion rates are relatively high (see Sec. 
5t and the protostar has an extremely large radius (e.g., 
Hosokawa et al.ll2012i ). In order to continue the calculations 
in these cases, we switch to an analytic stellar model, which 
agrees with the numerical results well, only when the accre¬ 
tion rate exceeds 1CU 2 Mq yr _1 (see also Appendix lAl). 


2.3 Classification of Star Formation Modes: 

Evaluation of Local Jlw 

To derive the primordial stellar mass distribution including 
both Pop III. 1 and III.2 d stars, we need to evaluate the local 
intensity of the photodissociating radiation at each cloud in 
a self-consistent manner with the spatial distribution of the 
primordial stars. For this purpose, we perform the following 
post-processing calculations: 

1. The primordial gas clouds are sorted by their formation 
redshifts z*i, the epoch when the cloud’s central density 
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Figure 3. Distributions of DM mini-haloes that host primordial star-forming clouds as functions of the formation redshift (left-hand 
panel) or the virial mass (right). In the left-hand panel, the bin width is A z = 1 and the different colours indicate the different ranges 
of the virial masses: M vir /M 0 < 10 5 ' 4 , 10 5 ' 4 < M vir /M 0 < 10 5 6 , 10 5 ' 6 < M vir /M 0 < 10 5 ' 8 , 10 5 ' 8 < M vir /M 0 < 10 6 °, and 
106.0 0 M v ; r /M 0 . In the right-hand panel, the bin width is logarithmically equal, AM v ; r (M) = (10 01 — 1) X M, and the different colours 
indicate the different ranges of the redshifts: 10 < z < 14, 14 < z < 18, 18 < z < 22, 22 < z < 26, and 26 < z < 30. The black lines show 
the sums over all virial masses or redshifts. The redshifts and virial masses are measured when riu C en = 10 7 cm -3 . 


reaches 10 7 cm -3 . We first begin with assigning the stellar 
mass to the halo which has the highest 2 * 1 . 

2. We then evaluate the local FUV intensity for the clouds 
forming at the lower redshifts. We assume that the star ra¬ 
diates at a constant luminosity Qlw during its lifetime t* 
and dies at redshi ft z*i = z * i — A z(t,) and adopt the nu¬ 
merical results of ISchaererl (120021 ). who present the values 
of Qlw an d t* as functions of the stellar mass. We evaluate 
the FUV intensity when a cloud is in the pre-collapse stage 
Zini (nH,cen = 10 cm -3 ). We count the FUV intensity from 
stars which have 2*1 > z- lrl \ > 2*2 to obtain the local field. 

3. The stellar mass in a halo forming at the lower red¬ 
shift is determined using the correlation derived from the 
2D RHD simulations. We classify the star as Pop III. 2d, 
if the local FUV intensity is higher than the critical value 
J 2 i,crit = 0.1. The mass of a star formed within the irradi¬ 
ated cloud is largely determined by the cloud’s properties 
for Pop III. 1 cases, whereas it is also dictated by the local 
FUV intensity for Pop III. 2d cases. 

4. We repeat steps 2 and 3 for the clouds one by one in the 
order of decreasing formation redshift. The stellar masses are 
accordingly determined by the above procedure. 

The intensity of the FUV radiation at the location of 
the ith cloud contributed by the star in the jth cloud is 

rj _ /esc hv avg Qlw..7 f A . 

LW - j “ 7T Ai/LW 47rd?. ’ W 

where / esc is the fraction of LW photons that escape from the 
jth halo, hv avg is the average photon energy emitted from a 
Pop III star in the LW band, A^lw is the difference between 
the maximum and minimum frequencies of the LW bandsQ 

4 favg = i'max = 13.6 eV (3.288 X 10 15 Hz) and ;z m i n = 11.2 eV 
(2.708 x 10 1S Hz). 



M star [M 0 ] 


Figure 4. Accretion histories of mass accreting protostars, which 
form from the same cloud (ID = 4 in Table [Til but with different 
FUV intensities J 21 = 0, 0.1, 1, and 10, as a function of the stellar 
mass. When the mass accretion rates fall below 10“ 4 M 0 yr“ 1 . 
we stop the simulations and assume the masses at this moment as 
the final stellar masses (50, 61, 883, and 496 M 0 , respectively). 


and dij is the distance betwe en the centres o f the ith and /th 
clouds jAgarwal et al.ll2012l ). According to iKitavama et all 
(l2004h . the escape fraction f esc should be around 0.1 — 1 for 
the mini-haloes considered with M v n = 2 x 10 s — 2 x 10® 
M 0 . We adopt the fiducial value of / esc = 0.5. The local 
FUV intensity at the ith cloud J 2 i,i is then calculated by 
summing up the contributions from all the nearby clouds 
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Figure 5. The local FUV intensity field, J 21 ? in the same comoving cosmological volume with (3 h ~ 1 Mpc) 3 at z = 25, 20, 19.5, 19, 
and 15. The colour contours indicate the FUV intensity ranging from J 21 = 0.025 — 6.3 (blue to red). The yellow and red clumps show 
the active Population III. 1 and III.2 d stars, respectively. The local FUV radiation field decreases with decreasing redshift (from left- 
to right-hand panels), because the typical stellar mass becomes lower and separations between the stars are stretched by the cosmic 
expansion. 


with 2 * 1 ,j Zixn,i ^- > Z* 2 ,j 5 
N 

JlW,! = jfw,i ■ (5) 

j 

We only consider the contemporary sources in the simulation 
box and ignore the background radiation (see Sec. 16.111 . 


3 OVERVIEW 

We present an overview of our main results in this section. 
Fig. [3](a) shows the redshift distribution of 1540 mini-haloes 
located in our cosmological simulation. Our sample ranges 
from z = 30 to 10 and peaks around z ~ 20. Fig. QJf b) shows 
the mass distribution of the haloes, which indicates a peak 
around M V - 1T ~3x 10 5 Mg. The samples in Paper I did not 
show such a peak, indicating some selection bias associated 
with the small cosmological volume and the limited statistics 
due to the small number of samples. In Paper I, we show that 
the halo mass is the key parameter to determine the stellar 
mass. Thus, the more accurate halo mass distribution of the 
current samples can be expected to yield a modified stellar 
mass distribution. 

Fig. [4]shows the mass accretion histories on to a proto¬ 
star, obtained from our 2D RHD simulations. As described 
in the previous sections, the computationally intensive RHD 
calculations are needed to examine the correlation between 
the final stellar mass and the intensity of the photodissoci- 
ating radiation. In all of the cases presented, the mass ac¬ 
cretion rate decreases sharply at some point: i.e., when the 
protostellar feedback shuts off mass accretion. In general, for 
the larger J 21 , the gas temperature in the accretion enve - 
lope is higher fe.g.. lOmukaifeoQll ; lO’Shea fc Normanll2008l) . 
causing higher accretion rates on to the protostar. Indeed, 
we find that the final stellar mass is larger for higher FUV 
intensities. We will account for this effect when estimating 
the final stellar masses for a given FUV intensity in order to 
construct the mass distribution of Pop III.2 stars (see Sec. 

EJ. 

Finally, we classify our data of 1540 primordial clouds 
into Pop III. 1 and III. 2d cases, according to our criteria 
based on the local photodissociating FUV intensity (J 21 < 
0.1 or not). Fig.[5]shows the spatial distributions of the FUV 



M star [MJ 


Figure 6. Mass distribution of Pop III. 1 and III.2 d stars. The 
dotted line shows the sum of the two populations (see also Fig. 
H7I below). 


intensity J 21 and the locations of active Pop III. 1 (yellow 
dots) and III. 2d stars (red dots) in the cosmological volume. 
There are small clusters of Pop III. 2d stars at high redshifts, 
where the photodissociating radiation is very strong. We de¬ 
termine the stellar masses from the cloud properties and the 
computed FUV intensity. Fig. E] shows the mass distributions 
of primordial stars including Pop III. 1 and III. 2d stars. The 
Pop III. 1 mass distribution has two peaks at ~ 250 and 
25 Mg. The two characteristic masses reflect different paths 
of thermal evolution during the early runaway collapse stage, 
which are driven by either H 2 cooling alone or H 2 cooling 
plus HD cooling (Paper I). Pop III. 2d stars have relatively 
large masses, clustering around ~ 400 Mg. We discuss the 
redshift dependence of the mass distribution in Section [ 6 ] 
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Figure 7. Scatter plots and histograms of basic properties of the 
gas clouds on the virial scale. The left-hand panels show, from top 
to bottom, spin parameters of DM, spin parameters of baryonic 
components, and the offset angles between angular momentum 
vectors of the two components. The right-hand panels show their 
histograms; the colours represent the same redshift ranges as in 
the right-hand panel of Fig. [3] 


4 COSMOLOGICAL SAMPLE OF 

PRIMORDIAL STAR-FORMING CLOUDS 

We study the physical properties of the 1540 primordial star¬ 
forming clouds found in our cosmological simulation. In this 
section, we ignore the effect of FUV fields and assume that 
all the clouds bear Pop III. 1 stars, allowing a direct compar¬ 
ison with the statistical results presented in Paper I. The de¬ 
pendences of cloud properties and stellar mass distribution 
on the formation redshifts are discussed. 

The cloud properties are calculated at two different 
mass scales by averaging over the gas in the virialized DM 
haloes and over the gravitationally collapsing gas at the 
Jeans scale. We define the boundary for the former to be 
the halo virial radius, within which the average matter den¬ 
sity is 200 times higher than the cosmic mean value. For the 
latter mass scale, the boundary is defined as the cloud ra¬ 
dius, where the ratio of the enclosed mass to the local Jeans 
mass (Eq. [l]) has its maximum value. We present the statis¬ 
tical analyses for the cloud properties calculated for the two 
mass scales. 



Radius [pc] 

Figure 8. Radial gas density profiles of primordial star-forming 
clouds when rtu C en = 10 7 cm -3 . The grey region shows the max¬ 
imum and minimum for all clouds at each radius. The coloured 
lines show averaged profiles for FD-cooling clouds at five red- 
shift ranges and the black line shows the averaged profile for 
HD-cooling cloud s. The dotted line sho ws the power-law slope, 
ng oc r -2 - 2 (e.g.. lOmukai fe Nishi|[l998l l. 


1 x 10 6 Mq. The virial temperature of a star-forming halo 
is about T v i r ~ 1000 K@, and thus the virial mass is 

/1 i - 3 / 2 

Afvirial.Sa ~ 4 x 10 s f-g-J Mo, (6) 

which gives M v i r = 2.1 x 10 5 and 9.9 x 10 s Mq at z = 
30 and 10, in good agreement with the typical masses of 
our mini-haloes. As shown in Fig. [3] the average halo mass 
increases with decreasing redshift, which is also consistent 
with the redshift dependence in Eq. ©. We have more than 
100 haloes per each redshift bin in the range of z = 22 — 14, 
which allows us to study even redshift dependences of the 
properties of our samples. 

One of the important quantities at the halo scale is the 
spin pa rameter, A = fvir/y2 R vir V v i r (following the defini¬ 
tion of iBullock et alj 200 il l , which characterizes the rota¬ 
tional degree of a halo. Fig. [7] shows the spin parameters for 
both DM and baryonic components and the relative angle of 
their angular momentum vectors. The log-normal distribu¬ 
tions and the time (redshift) evolution are consistent with 
previous studies. At high redshift, the baryon spin parameter 
is lower than that of DM. The distribution of baryon spin pa¬ 
rameter becomes close to that of DM at lower redshift (after 
z ~ 14) because of the momentum redis tribution between 
the two components (Ide Souza et al.|[2013l ). The average an¬ 
gle is Save — 35°, suggesting that the spin vectors of the two 
components are roughly aligned with each other in most 
haloes, although with a few exceptions whereby 6 > 90°. 


4.1 Virial Scale: DM Mini-Haloes 

The distributions of the formation redshift and the virial 
mass of mini-haloes are shown in Fig. [3] We see that most 
of the mini-haloes form at z = 30 — 10 with M v i r = 2 x 10 5 — 


5 This is the critical temperature for a primordial cloud to col¬ 
lapse with efficient H 2 molecular co oling which is weakly depen¬ 
dent on redshift ('e.g.. [Glove r| |2Q13T) . In fact, our samples of DM 
mini-haloes have temperatures close to the critical value (see also 
Paper I). 
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n H [cm -3 ] M(r) [M 0 ] 


Figure 9. Averaged profiles of primordial star-forming clouds when nn.can = 10 7 cm -3 in the same manner as in Fig. [8] (a) gas 
temperature, (b) H 2 fraction and ratio of HD to H 2 , (c) gas infall rate, and (d) Keplerian factor (/Kepler = D-ot /"^Kepler) • I 11 panel (b), 
the coloured lines show averaged H 2 fractions for H 2 -cooling cases at five redshift ranges. The dotted and dashed lines show averaged 
/hd//h 2 f° r ah HD-cooling and H 2 -cooling cases, respectively. In panel (c), the dotted lines show the fitting functions for the averaged 
infall rates at the protostellar formation and are characterized by {0.34, 0.34, 0.23, 0.17, 0.14} X (M e n C /M@) _ °' 7 Mq yr" 1 in order of 
higher- to lower-redshift groups. 


Interestingly, we find a trend that the offset angle and the 
spin parameter are anticorrelated, i.e., the gas and the DM 
components rotate differently in slowly rotating haloes. 


4.2 Jeans Scale: Gravitationally Unstable Clouds 

In a small halo, the collapsing gas is initially heated adi- 
abatically until H 2 molecular cooling overcomes the com¬ 
pression heating. When the H 2 fraction reaches the criti¬ 
cal value (/h 2 ,crit = a few xlO -4 ), the temperature begins 
to decrease rapidly with increasing density. We follow the 
collapse until the central hydrogen number density riH.cen 
reaches 10 7 cm -3 . Among our samples, the time duration 
of the cloud collapse, from riH.cen = 10 to 10 7 cnD 3 , in¬ 
creases with decreasing the formation redshift. For instance, 
the duration averaged over samples for five different redshift 
ranges are 29.5 Myr (10 < 2 < 14), 18.0 Myr (14 < 2 < 18), 
12.0 Myr (18 < 2 < 22), 7.9 Myr (22 < 2 < 26), and 5.5 Myr 
(26 < 2 < 30). The clou ds collapse more s lowly at lower red- 
shifts, as also found bv lGao et al.l (120071} . In the following, 
we study the evolution of clouds forming at different red- 
shifts. 

In the Pop III. 1 case, the thermal evolution proceeds in 
one of two different modes depending on whether or not HD 
molecular cooling is effective. In most cases, the evolution 
is driven by H 2 cooling alone, whereas HD cooling becomes 
efficient only under limited conditions. In a rapidly rotating 
cloud, for instance, H 2 cooling reduces the temperature to 
T < 200 K because of the slower collapse and hence less ef¬ 


ficient compressional heating. HD formation begins at these 
low temperatures and the resulting HD cooling further re¬ 
duces the temperature to the cosmic micr owave background 
(CMB) floor at Tcmb(z) — 2.73 (1 + z) K. iRipamontil (l2007fl 
study this mode of cloud collapse using ID hydrodynamic 
simulations and discuss its overall occurrence rate. We also 
find that this HD-cooling mode occurs in our 3D cosmolog¬ 
ical simulations (Paper 1)0 

We examine how efficiently HD cooling operates during 
the collapse in our samp les. Following previous studies (e.g., 
iMcGreer fe Brvanll2008l . Paper I), we define the HD-cooling 
case as occurring when the fraction ratio /hd//h 2 exceeds 
10 -3 at riH.cen < 10 6 cm -3 . Among our samples, we also find 
intermediate cases, whereby the onset of HD cooling is some¬ 
what delayed; i.e. /hd//h 2 increases sharply in the density 
range 10 6 cm~ 3 < riH.cen < 10 7 cm -3 and eventually exceeds 
10~ 3 . We classify our 1540 samples of primordial clouds as 
one of three possible cases, depending on /hd//h 2 during 
the collapse: H 2 -cooling cases (1186 cases), HD-cooling cases 
(151), and the intermediate cases (203), where the numbers 

6 Another possible route for HD- cooling primordial star- 
formation is suggested by som e studies l|Uehaj^r^Jjiutsuka|[200^; 
Shchekino^J^^^silie^J^OO^; iPrieto et al.l 120121 : lBovkioet'*al1 
20141: IPricto ct al.l l2014ll . In this scenario, the merging of 
multiple dark matter haloes induces the formation of shock- 
waves in which HD formation is enhanced. Such merging often 
yields large mass haloes wit h a io > 10 7 [(1 + z)/ 20] 2 Mq 
(IShchekinov fc Vasilievl 20061) , whereas our HD-cooling clouds are 
found in lower mass haloes with ~5 x 10 5 M© at z ~ 10 — 16. 
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in parentheses indicate the number of the corresponding 
samples. The H 2 -cooling cases show a broad range of for¬ 
mation redshifts: 10 < 0 < 14 (98), 14 < 2 < 18 (423), 
18 < 2 < 22 (474), 22 < z < 26 (168), and 26 < 2 < 30 
(23), which allows us to study the redshift-dependence of the 
physical properties of the primordial clouds. On the other 
hand, the HD cooling and intermediate cases appear only at 
low redshift (2 < 18). 

Fig. |H1 shows the averaged radial profiles of gas density 
when tiH.cen = 10 7 cm -3 . They are well fitt ed b y the power- 
law f unction p oc r -2 ' 2 in the envelope fe.g. lOmukai fe Nishil 
EE3 D The slopes of profiles look very similar, but, for any 
given radius, the clouds’ densities can differ by more than 
an order of magnitude (see the grey region). Such differences 
between the clouds will lead to differences of the gas infall 
rates on to a central protostar during the subsequent accre¬ 
tion stage. For the H 2 -cooling clouds, the densities tend to 
be lower with decreasing formation redshifts. The densities 
of the HD-cooling clouds tend to be even lower than for H 2 - 
cooling cases, a result of the lower temperature during the 
collapse as shown below. 

Fig. El displays the averaged thermal (panel a), chem¬ 
ical (panel b), and dynamical (panels c and d) properties 
of the collapsing clouds. Fig. Ela) shows the temperature 
distribution as a function of gas density, which reflects the 
variation of H 2 and HD abundances shown in Fig. EE(b). We 
find a small but systematic redshift-dependence of the ef¬ 
fectiveness of H 2 cooling: lower temperatures for lower red¬ 
shifts. This redshift-dependence comes from the fact that the 
collapse time-scale, which controls the thermal evolution of 
clouds, is actually longer for the lower redshifts on average. 
The collapse time-scale is set at an early stage of the cloud 
formation, and thus depends on the cloud’s physical prop¬ 
erties measured at the virial scale (see Sec. 14.41 later). For 
instance, the mean density at the virial scale is lower for the 
lower redshifts, which qualitatively explains the slower col¬ 
lapse at the lower redshifts (also see section 4.2.2 in Paper I 
for a more thorough discussion). For the HD-cooling cases, 
the higher /hd //h 2 leads to an even lower temperature. The 
typical mass of a gravitationally unstable cloud is given by 
the Jeans scale at the local minimum of the temperature pro¬ 
file, which is around several x 10 3 Mg for H 2 -cooling cases 
and several x 10 Mg for HD-cooling cases. 

Fig. EEc) shows the radial distributions of the instan¬ 
taneous gas infall rate Mcloud = 47rr 2 pi; ra d as a function of 
enclosed gas mass, extending well outside the Jeans mass. 
For the H 2 -cooling cases, the expected gas infall rates de¬ 
crease with decreasing formation redshifts. For the HD- 
cooling cases, the infall rates are about 10 times lower. We 
address this issue more in detail in the next subsection. The 
dotted lines show the expected profiles after further collapse. 
We extrapolate the profiles from the outer regions and get 
power-law profil es M c i ou d oc M° 7 which agree with previ¬ 
ous results (e.g., iGao et ai1l2007f ). 

Fig. EEd) displays the degree of the azimuthal rota¬ 
tion velocity v TO t normalized by the local Keplerian velocity 
UKepier = y/GM{r)/r, where M(r ) is the gas mass within a 
given radius r. We calculate v ro t by averaging the velocity 


7 The DM density profiles also represent power-law features 
which at r ~ 10 pc transition from hqm <x r~ 7 to r~ 15 . 
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Figure 10. Histogram of the gas infall rates measured at the 
Jeans scale when npj C en — 10 7 cm -3 . The lines with the different 
colours represent the different redshift ranges from high (red) to 
low redshift (blue). The top axis gives the stellar mass estimated 
from Mj eans by Eq. Q. 


perpendicular to the total angular momentum vector for the 
gas within r. All clouds have similar values clustered around 
/Kepler = 0.5, wh ich agr ees w ith the results of previous stud¬ 
ies ( e.g., Yoshida et al]|2006i '). 

iMcKee fe Tanl (120081 ') modelled the accretion histories 
as a function depending on two parameters: /Kepler and K , 
a measure of the entropy of the accreting ga s. We can com¬ 
pare their modelled accretion rates (fig. 9 in lMcKee fe Tail 
l2008lf and our results (Fig. [Ht) to determine the correspond¬ 
ing parameters. By assuming /Kepler = 0.5, the other pa¬ 
rameter K is estimated to be 1 ~ 2. They also estimate 
the final stellar masses when the photoevaporative mass- 
loss rate overcomes the accretion rate. The estimated final 
masses with above parameters range over 150 — 300 Mg, 
which are consistent with our results presented in later sec¬ 
tions. 


4.3 Distribution of Gas Infall Rates 

After the onset of collapse, the newly formed protostar grows 
in mass via accretion from the surrounding envelope. It is 
well known that the accretion history in this phase signifi¬ 
cantly affects the protostar’s evolution and the strength of 
the UV radiative feedback against the accretion flow. Since 
the stellar radiative feedback ultimately shuts off the mass 
accretion, the mean accretion rate is a key parameter that 
determines the final stellar mass. 

Fig. [TO] shows the histograms of gas infall rates. The 
black line represents the entire collection of samples cover¬ 
ing all formation redshifts. We find two peaks that corre¬ 
spond to the different thermal evolution during the collapse 
stages: the right peak at ~ 3 x 10 -3 Mg yr _1 is primar¬ 
ily associated with the H 2 -cooling case and the left peak 


© 2015 RAS, MNRAS OQO.[Tlf22l 







































Mass of Primordial Stars 11 


10" 4 10" 3 10" 2 
M vir [M Q yr 'j 



30 25 20 15 10 

Redshift 


Figure 11. Variation of the gas infall rate at the virial scale 
as a function of formation redshift and virial mass. The different 
colours depict the different mass infall rates, according to the 
colour scale at the top. The distributions of the filled circles and 
dotted lines represent the variation of the sample data and its 
fitting function Eq. ©• 


at ~ 10 -4 Mq yr -1 is primarily associated with the HD- 
cooling case. As expected from the redshift-dependence of 
the thermal evolution, the high-rate (high-mass) peak asso¬ 
ciated with H 2 -cooling cases is found at all epochs, whereas 
the low-rate peak associated with HD-cooling appears only 
at low redshifts. Thus, the infall rate distribution changes 
from a single peaked function at high redshifts to the bi- 
modal function for low redshifts, as shown by the coloured 
lines. 

In Paper I, with 110 cosmological samples of the primor¬ 
dial star-forming regions, we show that the large variation of 
the gas infall rates leads to a spread of final stellar masses 
ranging from 10 to 1600 Mq. The results of the present 
study, with 10 times more samples, show the existence of a 
characteristic value of the infall rate which depends on the 
formation redshift. Thus, we can expect that the stellar mass 
distribution should also have analogous redshift-dependent 
characteristic values. 


4.4 Estimating Pop III.l Stellar Masses 

Here, we revisit the correlations between the stellar mass 
Mm. i and gas infall rates, which have been examined in Pa¬ 
per I. We first begin with considering the infall rates mea¬ 
sured at the Jeans scale Mj eam . Paper I shows that the gas 
infall rates measured when riH.cen = 10 12 cm -3 can be used 
as an accurate estimator, via an empirical fit, of the final 
stellar masses. In the present study, in order to obtain a 
larger sample at a reasonable computation cost, we follow 
the cloud collapse only to riH.cen = 10 7 cm~ 3 . We there¬ 
fore need formulae to predict the Pop III.l stellar masses 


based on information available at riH.cen = 10' cnW 3 . To 
this end, we reanalyse the results of our simulations in Pa¬ 
per I and determine the correlation between the final stel¬ 
lar masses and the gas infall rates at the Jeans scale when 
u-H,cen = 10 7 cm" 3 . In general, the infall rates are lower 
than those at tiH.cen = 10 12 cm -3 , but they are still well 
correlated with the stellar masses fFig. 1121) . The new fitting 
formula is given by 


Mm = 250 M 0 


Afjeans 


2.8 x 10~ 3 M 0 yr" 


(7) 


The top axis in Fig. |TU] depicts the stellar mass calculated 
by this equation. In this diagram, the positions of the right 
and left peaks correspond to M* ~ 264 and ~ 24 Mq, which 
represent the typical masses of Pop III.l stars forming via 
H 2 -cooling and HD-cooling modes, respectively. 

Next, we discuss the dependence of the stellar mass on 
the halo properties. Paper I shows that the gas infall rate at 
the halo scale M v u is also correlated to the stellar mass, but 
the correlation is somewhat weaker. Our much larger sample 
size, however, enables us to derive more reliable correlations 
than in Paper I. We reconstruct the fitting formula at the 
halo scale, which depends on two parameters: the formation 
redshift and the halo mass (see Fig. 1111) . 


M v ir = 1.1 x 10 3 Mq yr 1 


1 + *YY M vir 

20 / V 4 x 105 M 0 



We interpret the redshift dependence of the gas infall rates 
(or the stellar masses) shown as the coloured lines in Fig. 1101 
as follows. For the H 2 -cooling cases, the infall rates at the 
halo scale are well correlated with those at the Jeans scale. 
On the other hand, such a correlation is not found for the 
HD-cooling cases, because the infall rates at the Jeans scale 
are within a narrow range of Mj ean s ~ 10 -4 Mq yr -1 . The 
relative scaling of M v i r and Mj ea ns are described by 

Mjeans = 3 x M vir for M v ir > 10 -3 M q yr -1 (H 2 ) , (9) 
= 10“ 4 M 0 yr" 1 < 10“ 3 M 0 yr -1 (HD) .(10) 


By substituting them into Eq. 0 , we obtain 

0.7 

Mm. i(h 2 ) = 264 Mq 


10 -3 M 0 yr - 


Afm.l(HD) — 24 Mq . 


( 11 ) 

( 12 ) 


Eq. dill) for Mni.i(H 2 ) is modified by eliminating M v i r with 
Eq. JHJ as 


A7m.i(H 2 ) — 282 M 0 


1 + z 
20 


M v 


4 x 10 5 Mq 


(13) 


Under the assumption that primordial star formation mostly 
occurs in 3<r mini-haloes with T v i r = 10 3 K, we evaluate the 
typical value of Mni.i(H 2 ) by substituting Eq. ([6j into (1131) . 


Mih.i(h 2 ) 


282 M q 


fl + z\ oeo5 

V 20 ) 


(14) 


This equation provides the typical stellar mass for the H 2 - 
cooling mode of 282 Mq at 1 + z = 20, which is consistent 
with our results shown in Fig. 1101 Eq. (1 141) is also consistent 
with the fact that, at low redshifts, the mass distribution 
peak moves to lower stellar masses. 
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ID 

RHD 

J 21 —0 
Estimate 

J21 

RHD 

=0.1 

Estimate 

J21 

RHD 

=0.316 

Estimate 

RHD 

T21 —1 

Estimate 

^21 = 10 

RHD Estimate 

1 

— 

12.5 

— 

872.6 

— 

352.3 

238.1 

368.3 


X 

2 

15.7 

23.7 

89.1 

93.6 

— 

286.4 

— 

449.8 

183.2 

203.7 

3 

25.7 

32.1 

— 

150.8 


X 

— 

809.5 

— 

2438.7 

4 

49.9 

49.4 

61.0 

100.7 

— 

425.7 

883.3 

433.7 

496.1 

373.6 

5 

— 

59.1 

152.0 

317.5 

177.7 

236.7 

126.6 

176.3 

— 

503.8 

6 

197.2 

142.2 

40.9 

33.0 

- 

608.2 

29.1 

70.6 


X 

7 

— 

328.2 

235.4 

268.5 

175.7 

544.4 

— 

522.7 

— 

473.8 

8 

— 

352.3 

60.1 

20.8 

— 

397.5 

— 

655.6 


X 

9 

291.1 

384.0 

323.7 

406.1 

— 

367.8 

— 

593.6 

618.8 

176.3 


Table 1. Column 1: Index of the clouds, Columns 2 — 11: Final stellar masses in the unit of Mq for different FUV intensities: J 21 = 
0, 0.1, 0.316, 1, and 10. The masses are determined by two different methods: the 2D RHD simulations (left) and from the estimating 
formula (Eq. [7] right). The underbars indicate that the cloud experiences significant HD cooling during the run-away collapse stage. 
The crosses indicate that cloud collapse is prevented by photodissociating radiation. The dashes indicate the cases without final stellar 
masses calculated by 2D RHD simulations. 



Mjeans [ M © > T 'l 


Figure 12. Final stellar masses obtained in 2D RHD simula¬ 
tions as a function of the gas infall rate at the Jeans scale when 
7®H,cen = 10 7 cm -3 . The open circles and squares represent 16 
cases with J 21 > 0 (Pop III. 2d stars) and 5 cases with J 21 = 0 
(Pop III. 1 stars), in which the 2D RHD simulations successfully 
follow the evolution until the final stellar masses are fixed. The 
crosses depict the reanalysed results from Paper I obtained for the 
same central density. The dotted line shows our fitting formula 
for Population III stars (Eq. [TJ. 


5 EVALUATION OF POPULATION III.2 D 

STELLAR MASSES 

In this section, we determine how the mass of a Pop III. 2d 
star is affected by external FUV irradiation by performing 45 
simulation sets. We first conduct 3D hydrodynamical simu¬ 
lations to follow the early cloud collapse under different FUV 
fields up to a central density riH.cen = 10 13 cm~ 3 . The pri¬ 
mordial cloud gravitationally collapses in 41 cases out of the 
total 45 cases. In the other 4 cases, gas condensation and the 
gravitational collapse are prevented because H 2 molecules 
are completely photodissociated (crosses in Table [lj. 

Then, we switch to 2D RHD simulations to study the 
mass accretion stage until the stellar UV feedback finally 


shuts off the mass accretion on to the protostar. The calcu¬ 
lations are finished in 21 cases of the 41 collapsed cases and 
the resultant stellar masses are shown in left-hand columns 
for each J 21 in Table [T] (HD-cooling clouds are marked by 
underbars). In the remaining 20 cases (dashes in Tabic [TJ, it 
was difficult to numerically follow the stellar evolution with 
time-dependent rapid mass accretion, so that their stellar 
masses are estimated by using the procedure described later 
in this section. 

Fig. ll2l shows that the final stellar masses obtained from 
the 2D RHD simulations are well correlated with Mj ea . ns . 
We find that the correlation for Pop III.2 d stars is similar 
to that for Pop III. 1 stars in Paper I (Eq. [TJ. This can be 
understood by noting that, for the Pop III. 2d cases, pho¬ 
todissociating radiation affects the thermal structure of the 
envelope only in the early stage of the collapse. Once the col¬ 
lapse proceeds and the density increases, gas self-shielding 
becomes effective and the subsequent evolution is unaffected 
by the external radiation. The evolution after the birth of the 
protostar is only modified with a different accretion history 
resulting from the different temperature in the envelope. 

Having derived the above correlation, we can estimate 
the stellar mass without computing the detailed long-term 
evolution after the birth of the protostar. Instead, we only 
follow the early collapse stage in 3D until the central density 
reaches 10' cm -3 . We determine the final stellar mass for 
the remaining 20 cases by using the calculated Mj eans and 
Eq. JTJ which are shown in the right-hand columns for each 
J 21 in Table [T] 

We see that, overall, external radiation increases the 
resulting final stellar masses, though the dependence can 
be quite complex in some cases. Table [T] shows that the 
stellar mass does not always increase with increasing FUV 
intensity. For the clouds with ID = 6 and 8, for instance, the 
stellar mass decreases somewhat from the III. 1 value, when 
a weak FUV intensity (J 21 = 0.1) is present. Interestingly, 
the thermal evolution during the collapse actually changes 
from an H 2 -cooling mode to the HD-cooling mode (marked 
with the overbars) for these cases: i.e. a weak FUV field 
sometimes triggers the HD-cooling mode. We shall examine 
the interesting cases in greater detail in the next sections. 
We first review the overall trend of our results, and then 
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Figure 13. Physical properties of collapsing clouds with different FUV intensity J 21 for the different cases ID = 4 (left), 8 (centre), 
and 9 (right) listed in Table [l] Line colours in each panel represent the different FUV intensity J 21 = 0, 0.1, 0.316, 1, and 10 from 
black to yellow. The top three rows show the evolution of the gas temperature (panels a), chemical abundances (panels b), and the ratio 
between the collapsing time and free-fall time /collapse — ^collapse/tff (panels c) measured at the cloud centre. The dotted lines in panel 
(a) show the Jeans scale as a function of the gas density and temperature. In the panels (b), the solid and dotted lines represent the 
number fraction of H 2 molecules /h 2 and the fractional abundance of HD molecules /hd //h 2 • The bottom three rows show the radial 
profiles of the gas infall rate (panels d), number density (panels e), and radial velocity (panels f) as functions of the enclosed mass for 
n H cen = 10 11 cm -3 . The dotted lines in panels (d) show the mass accretion histories seen in the 2D RHD simulations following the 
evolution after the birth of protostars. Note that the horizontal axis represents the protostellar mass for these cases. 

© 2015 RAS, MNRAS OOQ.mf22l 


















































14 Hirano et al. 





10 5 10 6 10 7 10 8 10 9 10 '° 

n H l cm ' 3 ] 


Figure 14. Density (colour and line contours) and velocity (arrows) distributions around the collapsing centre for the same cases as in 
Fig. 1131 The snapshots for njj cen = 10 11 cm 3 are presented. The box size is 0.1 pc on a side. The top and bottom panels show the 
slices through the X — Y and X — Z planes, whereby the Z-axis corresponds to the rotational axis. There is a blank panel for ID = 8 
with J 21 = 10, because the collapse is prevented by the photodissociating radiation in this case. 


examine the detailed evolution for the clouds ID = 4, 8, and 
9. Finally, we describe how these results are used to derive 
the mass distributions of Pop III. 1 and III. 2d stars. 

5.1 General Properties 

Fig. [J3]depicts the thermal and dynamical evolution during 
the cloud collapse. The columns represent ID = 4 (left), 8 
(middle), or 9 (right). For ID = 4 and 8, HD cooling be¬ 
comes efficient enough to alter the thermal evolution with 
the weak FUV intensity J 21 = 0.1. The decrease of the 
temperature at riH.cen < 10 7 cm -3 is explained by the in¬ 
crease of HD fraction and hence of the cooling rate (pan¬ 


els a and b). As described in Section 4.2, HD cooling be¬ 
comes efficient when the cloud collapse is slow. We calcu¬ 
late the elapsed time, tcoiiapse, between two snapshots with 
n-H,cen and 10 x riH.cen. Fig. 1131 c) shows that the time- 
scale ratio /collapse = tcoiiapse/fff exceeds 3 in the early 
stage@ We also see that the collapse is decelerated around 
n.H,cen ~ 10® cm~ 3 , when HD cooling becomes ineffective. 
In addition to the effect of HD cooling, the H 2 -cooling cases 


8 We have confirmed that the thermal evolution seen in the 3D 
simulations is well explained with the one-zone models using 
/collapse as a f ree parameter (see appendix C in Paper I). 
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show substantial variations of the temperature, depending 
on the FUV intensity (for instance, the ID = 9 case). We 
conclude that the accretion histories after the birth of a pro¬ 
tostar are also dependent on the local FUV radiation field. 

The bottom three rows in Fig. [13] show the dynam¬ 
ical properties of collapsing clouds at the point in time 
when riH.cen = 10 11 cm -3 . The instantaneous gas infall rate 
47rr 2 pu ra d (panel d) depends on the radial distributions of 
the density (panel e) and velocity (panel f). The variations 
seen in the profiles originate from the different thermal evo¬ 
lution during the collapse. Fig. md) also shows that the 
mass accretion rates observed in the 2D RHD simulations 
(dotted lines), which follow the subsequent evolution after 
the birth of a protostar (the abscissa is now the protostel- 
lar mass), follow the dependence on FUV intensity expected 
with the snapshots before the birth of the protostar (solid 
lines). 

Overall, the photodissociating radiation has a stronger 
impact on those clouds that evolve on low-temperature 
tracks for J 21 = 0 (see e.g. ID = 4 in Table [T] and Fig. 
EGO). In such cases, the thermal evolutionary path moves to¬ 
wards higher temperatures with increasing J 21 , which results 
in larger accretion rates and hence ultimately more massive 
stars. The effect is most significant for the clouds that would 
normally evolve on a HD-cooling path when FUV radiation 
is absent (e.g., ID = 1 — 5 cases). However, the dependence 
of the thermal evolution on the FUV intensity is not always 
simple nor even monotonic in some cases. For ID = 8, whose 
thermal evolution follows the typical H 2 -cooling mode for 
J 21 = 0, HD cooling dominates when J 21 is increased to 0.1. 
In the following, we look at individual cases in more detail 
to understand this complex behaviour. 

5.2 Individual Cases 

To study the evolution of individual cases, it is useful to ex¬ 
amine the 3D distributions of gas density and velocity for 
the different FUV intensities. Fig. [14] displays such distribu¬ 
tions in the central region with 1 pc on a side. We see that, 
overall, the morphology of a cloud changes from a disc-like 
shape to a more spherical one with increasing J 21 (from left 
to right). This can be understood as follows. At the higher 
FUV intensity, H 2 and HD molecules are photodissociated 
more efficiently, which makes the equation of state (EoS) of 
the primordial gas stiffen The collapse is decelerated in such 
a case and there is more time for angular momentum redis¬ 
tribution within the cloud. As a result, the cloud becomes 
more spherical. Once the density increases to the point that 
self-shielding effect becomes important, the collapse occurs 
rapidly with a weaker rotational support. This partly ex¬ 
plains why the strong FUV irradiation increases the gas in¬ 
fall rate. An exceptional case is for ID = 9 with J 21 = 10, 
which has significant rotation, resulting in a disc-like struc¬ 
ture with spiral arms. We describe the evolution of each case 
below. 


5.2.1 ID = 4 : disabling the HD-cooling path with FUV 
radiation 

We first consider the case ID = 4, for which the cloud evolves 
on the HD-cooling path for the Pop III. 1 case (J 21 = 0). 


Fig. [14] shows that the cloud rotates rapidly, which is also 
inferred from the low infall velocity shown in Fig. 113f fl. We 
remind the reader that the HD-cooling path emerges when 
the collapse is slow (Paper I), which could be the result 
of rapid rotation. The figures show that for J 21 =0.1 the 
evolution does not change much. The photodissociation of 
molecules hardly affects the evolution because of efficient 
self-shielding. The resulting stellar mass is almost the same 
as for the Pop III. 1 case. 

With the higher FUV intensities J 21 = 0.316 and 1, 
however, the thermal evolution tracks shift from the HD- 
cooling path to the H 2 -cooling path, because of the enhanced 
photodissociation of H 2 and HD molecules. The decreased 
coolants cannot cool the cloud sufficiently to follow a HD- 
cooling evolution. Once the thermal evolution begins to fol¬ 
low the H 2 -cooling path, the cloud becomes more spherical 
as described above. With the reduced rotational support, 
the collapse proceeds more rapidly and the gas infall rates 
are consequently higher than for the cases with weaker FUV 
irradiation J 21 < 0.1. 

The strong FUV irradiation with J 21 = 10 photodis- 
sociates H 2 molecules even up to 71 h = 10 9 cm -3 (see Fig. 
EM- With the reduced abundance of coolants, the cloud 
collapses slowly until three-body H 2 formation begins. The 
cloud has a spherical shape as a result of angular momentum 
redistribution over such a prolonged collapse time. The small 
rotational support results in the highest gas infall rates for 
M(r) < 100 Mq (or 71 h > 10 9 cm -3 ) among the five cases 
examined. In the outer part of the envelope, however, the 
infall rate is lower than for J 21 = 0.316 and 1, because of the 
slow collapse before self-shielding becomes effective. Due to 
the low infall rates in the outer low density regions (e.g., for 
riH.cen = 10 7 cnh 3 , see Fig. 1131 11. the final stellar mass with 
J 21 = 10 is lower than with weaker FUV irradiation. 

The ID = 4 cases represent the typical variation of 
the stellar mass with external photodissociating radiation: 
the stellar mass increases with FUV intensity J 21 . However, 
the trend saturates around J 21 ~ 10 and then the stellar 
mass decreases for even higher FUV intensity. Cloud collapse 
would ultimately be prevented with very strong FUV irra¬ 
diation, which is actually seen in some other cases (crosses 
in Table [[J. 

5.2.2 ID = 8: triggering the HD-cooling path with weak 
FUV radiation 

The next cases we examine in detail are for ID = 8, whereby 
for J 21 = 0 the cloud collapses along the H 2 -cooling path. In¬ 
terestingly, the HD-cooling path appears only for the weak 
FUV irradiation J 21 = 0.1 in this case. The main effect 
of the FUV radiation is photodissociating molecules. With¬ 
out molecules, the EoS becomes stiffer and the gas collapses 
slowly. However, the slow collapse actually promotes the for¬ 
mation of H 2 and HD molecules. Enhanced self-shielding 
eventually stops further photodissociation, triggering the 
HD-cooling mode in the subsequent evolution. At higher in¬ 
tensities J 21 = 0.316 and 1, the thermal evolution tracks 
return to the H 2 -cooling path, because the self-shielding is 
insufficient to prevent photodissociation. With J 21 = 10, the 
cloud does not collapse even after 10 8 yr. The strong pho¬ 
todissociation completely quenches star formation in this 
case. Our results suggest that the critical FUV intensity for 
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Figure 15. Stellar masses for the Pop III. 2 d cases with different FUV intensity J21. The left-hand panel shows the stellar mass of Pop 
III.2 d cases Mni.2 D (J21 > 0) plotted versus Pop III. 1 cases Mm .1 (J21 = 0), whose different values represent different gas clouds (for 
ID = 1 ~ 9 from left to right). The filled symbols indicate the stellar masses determined by 2D RHD simulations. The values marked 
with the open symbols, on the other hand, are calculated using the analytic formula (Eq. [7J The solid line represents the equal mass 
boundary Mm.l = Mni.2 D • The dotted lines show the fitting functions which follow the dependence of Mni.2 D on J21- The right-hand 
panel shows the correlation of Mni.2 D with J21. The symbols with error bars depict the averaged masses and variances of the same J21 
cases. The stellar masses for cases without 2D RHD simulations are given by the estimated values. The line shows the simplified relation 
Mm. 2 0 3 s a function of J21 (Eg. 1151 Afni.2 D = 100, 300, 500, and 400 Mg for J21 = 0.1, 0.316, 1, and 10), which is used for evaluating 
Pop III.2 d stellar masses. 


preventing collapse depends on the physical properties of 
the cloud such as its morphology and the degree of rotation. 
This is likely, because the strength of the self-shielding ef¬ 
fect depends on the radial profiles of the density within the 
cloud. 


5.2.3 ID = 9: reducing temperature with strong FUV 
radiation for the H 2 -cooling path 

Finally, we consider ID = 9 cases, which show an exceptional 
behaviour for strong FUV irradiation. The gas infall rate 
and resulting final stellar mass monotonically increase with 
increasing J 21 for J 21 < 1. For J 21 = 10, however, in spite 
of the intense FUV irradiation, the cloud evolves through 
the lowest temperature path among all cases examined (Fig. 
HSr). This is because the strong photodissociation deceler¬ 
ates the collapse even after the density becomes very high. 
In fact, Fig. 1131 c) shows that the time-scale ratio /collapse is 
around 5 at nu > 10 9 cnD 3 . Paper I finds that once three- 
body H 2 formation becomes effective, the gas temperature is 
lower for the longer collapse time-scales (see their appendix 
C); i.e., more H 2 molecules, the coolant, are formed with the 
slower collapse. The enhanced coolant abundance results in a 
softer EoS, which in turn lead s to a more flattened morpho l- 
ogy as the collapse advances (IHanawa &: Matsumotoll200d) . 
In fact, the cloud morphology accordingly changes from disc¬ 
like to a more spherical structure for the cases with J 21 ^ 1, 
but shows the highly rotating disc-like shape for J 21 = 10 
as seen in Fig. 1141 


5.3 Formula for Estimating Population III. 2d 
Stellar Masses 


We have seen that the masses of Pop III. 2d stars Mm. 2 D 
vary not only among different gas clouds but also by the 
effect of FUV irradiation. Despite the fact that Mm. 2 D is 
dependent on multiple parameters in a complicated way, 
the overall trend can be modelled by using the numerical 
simulation results obtained in this paper. Fig. 1151 a) summa¬ 
rizes how the stellar mass Mm. 2 D varies with the parameters 
{Mm.i and J 21 ). Note that, in the figure, the different Pop 
III. 1 stellar masses indicate that the gas clouds are differ¬ 
ent. Overall, the mass growth rate Mm. 2 D /Mm.i increases 
with decreasing Mm.i for a given FUV intensity; gas clouds 
with low temperatures are more susceptible to FUV irra¬ 
diation. On the other hand, the mass increase is relatively 
small for the cases with the higher Pop III. 1 masses, i.e., 
Mm.i > 100 Mq. In Fig. llSl a). we see that the same colours 
for a given J 21 are distributed around the similar values of 
Mm. 2 D . The Pop III. 2d stellar mass does not change radi¬ 
cally within a given FUV intensity J 21 . Table [I] shows the 
diversity of stellar masses for a given J 21 over a factor of 10 
but there is also a systematic dependence. Fig. 1151 b) shows 
the averaged values and variances for given values of J 21 . 
To obtain the averaged values and its deviations, we use the 
estimated stellar masses for the cases without the 2D RHD 
results. There is a general trend that Mm. 2 D increases for 
0.1 ^ J 21 ^ 1 but decreases gradually for J 21 ^ 1. We model 
such variations of Mm. 2 D for different J 21 by the solid line 
as 


Mm. 2 0 (J 21 ) 


'900 ■ 10 +0 - 9te (if 
500 • 10 +0 ' 44a: (if 
500 ■ 10~° lte (if 
.400 (if 


- 1 < x < --0.5), 


— 0.5 < x < 0), 
0 < x < 1), 
1 < x ), 


(15) 


© 2015 RAS, MNRAS OOO.ITH221 























Mass of Primordial Stars 17 


where x = log 10 (J 2 i). 


6 MASS DISTRIBUTION OF PRIMORDIAL 

STARS 

In this section, we calculate the mass distribution of Pop 
III. 1 and III. 2d stars using the non-biased cosmological sam¬ 
ples of the primordial star-forming regions (Sec. [4]). We de¬ 
termine the stellar mass for each halo by considering the 
physical properties of the cloud for Pop III. 1 cases (Eq. [71) 
and intensity of the local FUV radiation field for Pop III. 2d 
cases lEq.fTol). 

6.1 Classification into Population III.l and III.2 d 
S tars 

In the cosmological context, as shown in Fig. [0 the ear¬ 
liest generation of stars are formed in filaments or knots 
of the large scale structure and hence the star-forming re¬ 
gions are distributed in a biased manner. Radiation emitted 
by stars formed early affects subsequent star formation in 
nearby regions. As described in Section 12.31 we evaluate the 
local FUV intensity at each cloud considering the spatial 
distribution of other stars at the moment when the cloud 
central density reaches 10 cm” 3 . We assume that Pop III. 2d 
stars form, when J 21 > 0.1. Fig. [7]shows the resultant FUV 
radiation fields in the cosmic volume at five different red- 
shifts. The mean FUV intensity decreases with decreasing 
redshift mostly because spatial separations between stars are 
stretched by cosmic expansion. At z = 25, there is a small 
cluster of Pop III. 2d stars that are formed under the influ¬ 
ence of a strong FUV field. Once the bright massive stars 
form in a dense region, the resulting FUV field could trig¬ 
ger the sequential formation of Pop III. 2d stars in nearby 
haloes. At z = 15, on the other hand, there are no Pop 
III. 2d stars in the simulation volume because of the weak 
FUV radiation field resulting from the large physical sepa¬ 
ration between star-forming regions. 

Fig. [16] shows the histogram of the calculated FUV 
intensity and its redshift dependence. About one third of 
haloes is irradiated by a strong FUV field and thus satisfies 
our criterion for Pop III. 2d star formation. At z > 20, the 
clouds have a large variation of J 21 , ranging from 0.01 to 
10; about half of the clouds are exposed to FUV radiation 
with J 21 > 0.1. At z < 20, however, the fraction of III. 2d 
stars rapidly decreases. The transition is caused by a com¬ 
bination of the increase of physical separations between the 
clouds, the decrease of the number of FUV-active stars, and 
the systematic decrease of stellar masses by the promoted 
HD-cooling mode of star-formation. We also find that, at 
15 < z < 25, almost all the clouds have non-zero intensity 
J 21 > 0.01. However, such a weak FUV irradiation causes 
negligible effects on the cloud collapse. 

Regarding clouds under strong FUV irradiation with 
J 21 > 10, we have 5 such samples at z = 21 — 26 with in¬ 
tensities of J 21 = 12.1, 14.1, 17.7, 70.3, and 97.3. However, 
these are still below the critical value of J|i' 4 oc 0(1000), 
above which the so-called direct collapse mig ht be triggered 
in at omic-cooling haloes (e.g., IWolcott-Green fe Haimanl 
l201lh . In the current cases with the weaker FUV radiation 
and lower mass haloes, the gas clouds likely yield massive 



Figure 16. Histogram of the normalized intensity of photodis- 
sociating radiation J 21 at the position of each primordial cloud 
when njj C en = 10 cm” 3 . The different colours depict the same 
redshift ranges as in the right-hand panel of Fig. [3] The dotted 
vertical line represents the critical value of J 21 = 0.1, above which 
Pop III.2d stars form. The haloes with J 21 > 10 and J 21 < 10” 3 
are included in the rightmost and leftmost bins, respectively. 


Pop III.2 d stars but not the massive BHs as prospected in 
the direct collapse scenario. 

In the next section, we determine the stellar mass for 
each of 1540 star-forming clouds according to the star- 
formation mode: we use Eq. 0 for Pop III.l stars and Eq. 
(nsj for Pop III.2 d stars. 

6.2 Mass Distribution of Primordial Stars 

Fig. [TT] displays the mass distributions of Pop III.l and 
III. 2d stars. The black solid lines show the mass distribu¬ 
tions integrated over redshifts. As expected from Fig. 1101 the 
Pop III.l mass distribution has two peaks around M* ~ 250 
and 25 Mg, which reflect the contributions of the H 2 -cooling 
and the HD-cooling modes (see Sec. 14.311 . The mass dis¬ 
tribution of Pop III. 2d stars is shifted to larger masses 
M* ~ 400 Mq. We find a wide mass range for each popu¬ 
lation: 50 < Mm.iH 2 /M 0 < 1000, 10 < tfm.i HD /M 0 < 50, 
and 100 < Mm. 2 D /M 0 < 1000, respectively. About a half of 
M* > 200 M 0 stars are Pop III. 2d stars. The mass distribu¬ 
tions for the different redshifts are represented by lines with 
different colours. For Pop III.l cases, as suggested by Eq. 
da, the high-mass peak, corresponding to the H 2 -cooling 
mode, gradually shifts to lower stellar masses with decreas¬ 
ing redshift; from 375 M 0 at z = 30 to 191 M 0 at z = 10. 
The mass distribution of Pop III. 2d stars shows a similar 
trend, which is caused by the decrease of the average FUV 
intensity as shown in Fig. [16] We find no such shift for the 
HD-cooling cases, probably because of the small sample size 
at lower redshifts. Overall, the decrease of the Pop III. 2d 
fraction, the decrease of Pop III.1 h 2 stellar masses, and the 
increase of low-mass Pop III.Ihd fraction explain the shift 
of the distribution to the lower masses with decreasing red¬ 
shift. 

Because we use the analytic function Eq. 0 to estimate 
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Figure 17. Resultant mass distributions of Pop III. 1 (left) and III.2 d (right) stars for the different redshifts. The different colours 
represent the same redshift ranges as in the right-hand panel of Fig. [3] The black solid lines show the total distributions over all redshifts 
for each population whereas the dotted lines show the sum of them. 


the Pop III. 1 stellar masses from the gas infall rates at the 
Jeans scale, the mass distribution reflects that of the infall 
rates. Interestingly, except for Pop III. 1 cases formed via 
the HD-cooling mode, the stellar mass distributions are well 
described by power-law functions at both the low-mass and 
high-mass ends. For example, the mass distribution of Pop 
III.1 h 2 stars is approximately proportional to M*' 5 at the 
low-mass end and to M* at the high-mass end. This allows 
us to define the normalized stellar mass function, 


®(M») 


1.62 

M p /M 0 

1.62 

Mp /Mg 


/ m \ 2,5 

(_ij ,o r «,<«„, 

(nr for m ->• 


(16) 

(17) 


where M p is the peak mass given by Eq. m and 
is normalized by f >1 '(M*) dM * = 1. This is obviously quite 
different from the well-known Salpeter function. The above 
equation allows us to model the time-dependent mass func¬ 
tion of primordial stars. 


6.3 Star Formation Rate Density 

Fig. 1181 a) shows the star formation rate densities (SFRD) 
as a function of redshift. The primordial SFRD rises until 
2 ~ 20 and decreases afterwards (as Pop II star formation 
becomes the main mode). This evoluti on of the total SFRD 
is consistent w ith previous studies (e.g. jAgarwal et alJboij : 

1 Johnson et al.l [20131 ). but our results show clearly the con¬ 
tribution of Pop III. 2d stars. Significant Pop III. 2d star for¬ 
mation occurs after Pop III. 1 stars are formed and emit co¬ 
pious amounts of FUV photons. Remarkably, SFRDih.2 d 
approaches the same level as SFRDm.i at z > 20. At 

2 < 20, the fraction of Pop III. 2d stars rapidly decreases 
as shown in Fig. 1161 Note, however, we may be underes¬ 
timating SFRDiii.2 d at this point, because we ignore the 



Figure 18. Redshift evolution of SFRD of Population III. 1 and 
III.2 d stars (panel a) and the averaged J 21 for clouds with its 
variance (panel b). In panel (a), the dotted line represents the 
total of them. In panel (b), the grey dots show the scatter of J 21 
at each cloud and the dotted line represents the critical value of 
J 21 = 0.1, above which Pop III.2 d stars form. 


contributions of Pop II stars to the local and global FUV ra¬ 
diation field. The FUV background radiation intensity may 
exceed J 21 = 0.1 at z < 10 and the P op II SFRD may 
dominate at z < 15 (e.g., lAgarwal et al.l[20121) . If a large 
number of Pop II stars are formed during the epoch we con- 
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Fate 

Mass Range 
(Mq) 

III.l 

(10 < 

III.2d 
z < 14) 

(14 < z 

< 18) 

(18 < 

z < 22) 

(22 < 

z < 26) 

(26 < 

z < 30) 

Total 

NS 

M* < 25 

2.4% 


0.7% 








3 . 1 % 

LMBH 

25 < M* < 80 

4.2% 


2.3% 


0.3% 






6.8% 

PPISN 

80 < M* < 120 

1.4% 


2.7% 


0.7% 

0.1% 

0.1% 




5.0% 

PISN 

120 < M* < 240 

3.7% 


12.1% 

1.4% 

4.8% 

2.2% 

0.7% 

0.1% 



25.0% 

HMBH 

240 < M* 

2.0% 

0.1% 

13.8% 

4.0% 

13.0% 

14.2% 

5.6% 

5.7% 

0.9% 

0.8% 

60.1% 

Total 


13.7% 

0.1% 

31.6% 

5.4% 

18.8% 

16.5% 

6.4% 

5.8% 

0.9% 

0.8% 

100.0% 


Table 2. Column 1: Final fate of stellar evolution: neutron star (NS), low-mass black hole (LMBH), pulsation pair-instability supernovae 
(PPISN), pair-instabilit y supernova e (PISN), and high-mass black hole (HMBH), Column 2: Corresponding mass ranges at ZAMS for 
each fate from fig. 12 in I Yoon et alj ll20l2 i for the non-rotating case, Columns 3-12: Fractions of each fate in the sample for each redshift 
range, Column 13: Sum over each row. 


sider here, the number fraction of Pop III.2 d stars could be 
enhanced and the mass distribution at z < 15 may be sig¬ 
nificantly modified. There are a number of M, < 100 Mq 
Pop III. 1 stars forming via the HD-cooling mode, which is 
easily changed to the Fh-cooling mode even for even weak 
FUV fields (see Sec. [5J. The enhancement of the stellar mass 
by photodissociating molecules is significant for these cases 
(Fig. 1151) . However, because the total primordial SFRD de¬ 
creases for £ < 20, the uncertainties described above would 
not greatly change the overall mass distribution integrated 
over redshifts. 

Fig. I18l b) shows the redshift evolution of averaged local 
FUV intensity at each cloud. The averaged value decreases 
with decreasing redshift and falls below the critical value 
at z ~ 20, which is consistent with the decline of SFRD 
for Pop III.2 d cases (Fig. I18h ). I n comparison to the back¬ 
ground FUV f i eld ca lculated by lAgarwal et al.l (l2012fl and 
I Johnson et alj (|2013h . the local J 21 is at the same levels for 
s > 15 but starts to decline earlier for z < 15. This ear¬ 
lier decline is because of our ignorance of Pop II stars for 
calculating the FUV radiation fields. 

7 DISCUSSION 

7.1 The Final Fate of Primordial Stars 

The derived stellar mass distri bution can be use d to pre¬ 
dict the fates of primordial stars. lYoon et all (|2012|) perform 
stellar evolution calculations and categorize the final fates 
of primordial stars as a function of the stellar masses. We 
classify our data of primordial stellar masses into five dif¬ 
ferent populations depending on their final fates for non¬ 
rotating cases: neutron stars (NSs; M*/Mq < 25), low- 
mass black holes (LMBHs; 25 < M*/Mq < 80), pulsational 
pair-instability supernovae (PPISN; 80 < M,/Mq < 120), 
pair-instability supernovae (PISN; 120 < M*/Mq < 240), 
and high-mass black holes (HMBHs; 240 < M*/Mq). Ta¬ 
ble [2] summarizes the respective number fractions and their 
redshift-dependence. Because the mass distribution itself 
evolves over time, the relative occurrence of the various fates 
of primordial stars also changes over cosmic time. 

First of all, a large fraction of the primordial stars 
end their lives as HMBHs. They do not contribute to the 
early chemical evolution but leave remnant black holes that 
might have seeded the formation of supermassive black 


holes dLi et al.ll2007l ). A fraction of Pop III.1h 2 stars die as 
PISN that chemically pollute the surrounding inter-stellar 
medium, leaving the peculiar elemental abundance patterns 
associated with this type of supernovae. It has been a long¬ 
standing puzzle, however, that such distinct abundance pat¬ 
terns are not obs erved in Galactic metal-poor stars. Very 
recently, however, lAoki et~al1 d2014l ) report the discovery of 
a metal-poor star whose elemental abundance patterns can 
be explained by a supernova explosion of a very massive star 
with M, > 100 Mq. Intriguingly, they also estimate from 
their data obtained from the sloan extension for galactic un¬ 
derstanding and exploration (SEGUE) survey that such pe¬ 
culiar stars are rare, comprising only a few percent of metal- 
poor stars, in reasonable agreement with our simulation re¬ 
sult. One explanation for the challenges in finding PISN sig¬ 
natures may be that mini-haloes hosting such massive pri¬ 
mordial stars have masses too low to retai n the halo gas af¬ 
ter a n energetic supernova explosion (e.g.. lCooke fe Madaul 
120141) . Heavy elements produced in the stellar interior might 
be expelled to intergalactic space and are therefore absent 
from the gas out of which the metal-poor stars form. At 
lower redshifts, the Pop III.Ihd stars with M* ~ 25 Mq 
die leaving LMBHs or NSs. Nucle osynthesis in core-collaps e 
supernovae or hypernovae (e.g., lUmeda fc Nomotol I2003I ). 
which are expected for such progenitors, are likely to domi¬ 
nate in the early Galactic chemical evolution. 

7.2 Uncertainty in Mass Estimation 

In the present investigation, we have used simple analytic 
formulae for estimating the stellar mass from local properties 
such as the gas infall rate Mj eans and the FUV intensity J 21 . 
It would be important to improve the estimates by including 
physical effects that have not been considered here. 

For instance, recent 3D numerical simulations show that 
a circumstellar disc becomes grav i tationally unstab l e and 
fragments (e.g., Clark et al.l 1 20 111 : iGreif et all 1201 ll . 120121 : 
IStacv fc Bromrr] 20131 ). This could lead to the formation of 
multiple stellar systems rather than a single star in each 
mini-halo. In this case, the available gas in the envelope is 
divided among multiple protostars, so that the stars in such 
a multiple system could have relatively lower masses than in 
the c ase of a single star (e.g.. IPeters et aj]|2010l : ISusa et all 
120141) . However, there is an opposite effect; a large frac¬ 
tion of the protostars could rapidly migrate inward due to 
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clouds that were previously ionized by radiation from nearby 
stars. In principle, one can simulate the evolution of such 
clouds and nearby Hu regions by s olving the transfer o f 
EUV (hu > 13.6 eV) radiation fe.g.. lYoshida et al.ll2007h . 
This will likely add an additional p opulation of lower mas s 
primordial stars (e.g., M* = 17 Mm: lHosokawa et al.ll201~2i l. 
which can be verified in future studies. However, for the 
stellar mass distribution at the redshifts of interest in this 
paper (10 < z < 30), Pop III.2i stars are a minor component 
and do not affect the conclusions of this study. 


Figure 19. Same as Fig. [6] but calculated by decreasing the 
estimated stellar masses obtained from Eq. 0 by factors of 0.5 
(panel a) and 0.2 (b). 


gravitational torque, res ulting in frequent stellar mergers 
at th e cloud centre fe.g.. iGreif et al.ll2012l : IVorobvov et al] 
l2013ll . This more efficient accretion could enhance the for¬ 
mation of massive stars. The final stellar masses are deter¬ 
mined by the balance of the above competing effects, which 
would also be affected by s tellar radiative feedback (e.g., 
IStacy et alj l20lilSusall20ll). Moreover, recent magnetohy¬ 
drodynamic simulations show that even a small amount of 
turbu lence can ampli f y magnetic fields by a dynamo mecha- 
nism (ISur et al]|20ld : iFederrath et aLlfeoill ; ISur et akll2012l : 

I Turk et al.l 20121 1 . The existence of magnetic fields would in¬ 
crease angular momentum transfer in the disc via magnetic 
braking, preventing disc fragmentation as w ell as enhancing 
accre tion rates on to the protostars (e.g., iMachida fc Doil 
l2013l l. Simulating the long-term evolution including all of 
these effects is still challenging, but should be tackled in 
future studies. 

The number fraction of Pop III. 2d stars is dependent 
on the estimated mass of the Pop III. 1 stars, which in turn 
determines the strength of the FUV radiation field. It is pos¬ 
sible to see how critical this effect is by artificially modifying 
the scaling relation Eq. ©. Fig. ll9l shows the mass distribu¬ 
tions obtained using modified scaling relations, whereby the 
masses obtained from Eq. 0 are multiplied by either 0.5 or 
0.2. The local FUV intensity decreases because of the lower 
stellar masses, and the fraction of the clouds with J 21 >0.1 
that form Pop III. 2d stars also decreases. In panel (a), the 
numbers of the two populations are still comparable at the 
high-mass end, whereas in panel (b), the overall mass distri¬ 
bution is dominated by Pop III. 1 stars. Uncertainties in the 
masses of the Pop III. 1 stars affect the masses of the Pop 
III.2 d stars in this manner. 

There are also uncertainties regarding the local FUV 
intensity. In this paper, we only have counted photodissoci- 
ating radiation emitted by primordial stars. However, there 
could be other FUV emitters such as Pop II stars. Pop II 
stars dominate the cosmic star formation rate at low red- 
shifts. As noted in Section [2] however, including their FUV 
radiation should not largely modify the stellar mass distri¬ 
butions because a large fraction of the primordial stars form 
at sufficiently high redshifts. 

There is another type of Pop III.2 star that is formed 
under the influence of external radiation, the so-called Pop 
III.2i star (see also Sec. 0. We have indeed found potential 
sites of Pop III.2i star formation, i.e., stars forming in gas 


7.3 Cosmological Samples 

There is a variety of physical processes involved in the for¬ 
mation of primordial stars. Most of the uncertainties we have 
discussed are attributed to the complexity of the star for¬ 
mation process itself, in particular, during the stage of mass 
accretion on to the embryo protostar and the protostar’s UV 
feedback on the collapsing parent cloud. By contrast, the 
early evolutionary phases until the formatio n of a primor¬ 
dial gas cloud is relati vely well understood (lYoshida et al.l 
120081 : IGreif et aLll2012f ) and our results provide statistics of 
the physical properties of the cosmological primordial gas 
clouds, such as the radial profiles of density, temperature, 
and velocities. In future, we shall discuss the variety of cloud 
properties and resultant stellar mass. Furthermore, we find 
that there is a certain normal initial condition of the pri¬ 
mordial star formation and its redshift-dependence. 

Previous studies adopted different initial conditions and 
used different numerical codes with different physical pro¬ 
cesses implemented. It is thus often difficult to identify 
the origin of differences between studies. Direct compari¬ 
son would be easier if the same initial conditions are used. 
Our sample of haloes can provide a convenient test set as 
well as statistical properties of haloes under normal initial 
conditions [f] 

Ultimately, we will be able to determine the final stel¬ 
lar masses more accurately, when we conduct detailed 3D 
radiation magnetohydrodynamic simulations, including the 
additional physics described in Section 17.21 Alternatively, 
we can do the same more efficiently with knowledge of the 
distribution of initial conditions. Instead of repeating the 
entire multi-step procedure performed in the present study, 
we only need to update the estimation formula of the final 
stellar masses such as Eq. 0 by simulating the evolution for 
several representative cases. In this manner, the primordial 
stellar mass distribution can be easily reconstructed from a 
very large cosmological simulation. 
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APPENDIX A: ANALYTIC MODEL FOR THE 
ACCRETING SUPERGIANT PROTOSTAR 


The evolution of mass a ccreting protostars gre atly depends 
on the accretion rates. IHosokawa et al.l (12012} 1 show that, 
with high accretion rates of M > 0.04 Mq yr -1 , the pro- 
tostellar evolution qualitatively differs from that at lower 
rates. At these high accretion rates, the stellar radius mono- 
tonically increases with increasing the stellar mass, exceed¬ 
ing 10 3 Rq for M* > 100 Mq. The UV radiative feedback 
from such “supergiant protostars” is weak, because the re¬ 
sulting effective temperature is less than 10 4 K. In the cur¬ 
rent work, we frequently see this evolutionary stage in our 
2D RHD simulations of Pop III.2 d cases, where the accretion 
rates are relatively high. However, it is sometimes difficult 
to construct stellar models of these strongly accreting super¬ 
giant protostars by numerically solving the stellar structure 
equations, especially with highly variable mass accretion his¬ 
tories. 

Instead, we shall consider an analytic model for the evo¬ 
lution of the supergiant protostar. It is known that, during 
this evolutionary stage, the evolut ion of the stellar radius is 
well described by equation (11) in IHosokawa et ahl (120121 '). 

/ M* \ 1/2 

_R*(M») ~ 260 R 0 , (Al) 

a formula derived analytically assuming the Eddington lu¬ 
minosity and a stella r surface temperature T e g = 5000 K. 
As shown in fig. 5 of IHosokawa et ahl ( 2012 ?). however, the 
stellar radius obtained using a detailed stellar evolution code 
is slightly less than that predicted by Eq. ED- Neverthe¬ 
less, the functional dependence with an exponent 1/2 is well 
matched. We therefore retain the functional form of this de¬ 
scription but multiply the resulting radius by 0.8 to better fit 
the numerical results in Paper I, which also used a detailed 
stellar evolution code. The stellar effective temperature is 
nearly constant at T e g = 5000 K during this phase (see fig. 
12 in Paper I), and the stellar luminosity is calculated from 

L,(R*,T eff = 5000 K) = AirRlaTeg , (A2) 


where a is the Stefan-Boltzmann constant. As long as the 
mass accretion rate is above the critical value 0.04 Mq yr -1 , 
we calculate the stellar radius and luminosity using the 
above equations. Once the accretion rate falls below this 
value, we switch to another analytic model of the “oscillat¬ 
ing protostar” used in Paper I. 

After the mass accretion rate falls below 4 x 
10 -3 Mq yr -1 , the star begins to contract on a Kel vin- 
Helmholtz (KH) time-scale fe.g.. lOmukai fc Pallall2003h . 


Ikh 


Gtf, 

R*L* 


(A3) 


The changes of stellar radius and luminosity during a 


timestep dt are written as 


new 

-it* 

new 


R* + -—- • [R(M*) zams — R*] , 

CKH 

(A4) 

L* + -—- • [L(M*)zams — -£/*] , 

CKH 

(A5) 


where the properties of zero-age main sequence (ZAMS) 
stars are 

/M \°' 58 

R(M t ) zams = 3.109 x 10 _1 , (A6) 

/ TUT \ 130 

L(M,)ZAMS = 3.939 x 10 3 . (A7) 

If Ea.lAflvields a value R* new (M.) < Rzams , we use another 
relationship to get the stellar properties at the next step 


new 

new 


R* + - — -R(Mt) ZAMS 

IKH 



L, 


+ - - L{M„) ZAMS 

IKH 



(A8) 

(A9) 


By using the above-mentioned formulae for protostars 
going through the supergiant phase, we can compute the 
mass growth of a rapidly accreting Pop HI.2 d star. 
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